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A theoretical interpretation of the recent experimental studies of temperature evolu- 
tion in the course of time in the freely-expanding ultracold plasma bunches, released from 
a magneto-optical trap, is discussed. The most interesting result is finding the asymptotics 
of the form Tg oc t~(^-^^°-^) instead of which was expected for the rarefied monatomic 
gas during inertial expansion. As follows from our consideration, the substantially de- 
celerated decay of the temperature can be well explained by the specific features of the 
equation of state for the ultracold plasmas with strong Coulomb's coupling, whereas a 
heat release due to inelastic processes (in particular, three-body recombination) does not 
play an appreciable role in the first approximation. This conclusion is confirmed both by 
approximate analytical estimates, based on the model of "virialization" of the charged- 
particle energies, and by the results of ah initio numerical simulation. Moreover, the 
simulation shows that the above-mentioned law of temperature evolution is approached 
very quickly — when the virial criterion is satisfied only within a factor on the order of 
unity. 

PACS: 52.25.Kn, 52.27.Gr, 52.65.Yy 
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1 Introduction 



Apart from the previously- known types of nonideal plasmas [1], an active study of one 
more kind of the nonideal Coulomb's systems — the bunches of very rarefied ultracold 
plasmas created by laser cooling and ionization in the magneto-optical traps — was started 
in the recent years {e.g., review [2]). They are the classical (non-quantum) gaseous systems 
with a characteristic temperature of a fraction to a few Kelvin, whose Coulomb's coupling 
parameter T = e'^n^^^ /ksT can reach considerable values. For example, the immediately 
measured values for ions Fj are 2-^3 [3]; the corresponding estimates of Fg for electrons 
are less accurate and model-dependent, but they also give the values comparable to unity. 

The possibility of existence of such metastable plasma states was theoretically pre- 
dicted many years ago (see, for example, article [4] and references therein), but they were 
created experimentally only after a sufficient development of the laser cooling technique [5, 
6]. In the most recent time, similar systems began to be studied also by gas-dynamic cryo- 
genic installations [7]. Besides, creation of the same plasma states by artificial release of 
gaseous clouds from spacecraft was discussed long time ago [8, 9]. Unfortunately, the 
diagnostic possibilities in space still remain too limited to give reliable conclusions about 
the properties of the resulting plasmas. 

One of the most interesting results of the laboratory experiments performed by 
now was studying a temporal behavior of the temperature in the ultracold plasma clouds 
released from a magneto-optical trap and expanding freely in space. It was found, firstly, 
that the law of decay of the electron temperature at large times becomes "universal", i.e. 
independent of the initial conditions. (For example, even when the initial temperatures 
were scattered by 30 times, the values of Tg after a few microseconds deviate from each 
other by only 2-^3 times, and even less later [10].) Secondly, which is more interesting, the 
measured asymptotics had the form Tg oc t-(i-2±o.i) ^ ^^-i j^j^j^j jng^ga^j of which should 
be expected for the ideal rarefied gas without the internal degrees of freedom (7 = 5/3) at 
the inertial stage of expansion {i.e. when the plasma cloud expands with a constant rate, 
so that its size increases linearly in time, R{t) oc t). 

The most evident way to explain the substantially decelerated decay in the electron 
temperature is to take into account a heat release due to recombination of the charged 
particles. In the particular case of atomic ions, the most efficient channel is the three-body 
process, A'^ + e + e-^A + e, when one electron is captured by the ion, and the second 
electron carries away the excessive energy. Unfortunately, the recent attempts of quanti- 
tative modeling the observed law of temperature variation due to the heat release by the 
tree-body recombination were unsuccessful: the resulting dependence Te(t) differed only 
slightly from the adiabatic case (see, for example, the inset to Fig.3a in paper [lljl). 



^ Let us mention that the method of drawing the plots in paper [11] is somewhat confusing. According 
to the physical sense of the problem, the various laws of evolution should be confronted with each other 
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The aim of the present article is to show that all the experimentally measurable 
features in the temperature evolution (namely, both the appearance of asymptotics almost 
independent of the initial conditions and its particular form, close tot~^) can be naturally 
explained by the model of "virialization" of the charged-particle energies, i.e. actually 
by changing the equation of state of ultracold plasma for the case of strong interparticle 
Coulomb's interaction. As a result, it becomes unnecessary in the first approximation 
to take into account any inelastic processes, such as the heat release by the three-body 
recombinationFI 



We shall present below the estimates of the law of temperature evolution during expansion 
of cold nonideal plasmas that were actually performed over 10 years ago, before the 
experimental measurement of temperature in the magneto-optical traps [8, 9]. These 
estimates were done for some kinds of the plasma outbursts important in astrophysical 
applications. We are not going to discuss here these applications but would like to remind 
the basic calculations and the results obtained. 

First of all, since a kinetic energy of thermal motion of the monatomic ideal gas 
during its inertial expansion decreases as t^^, while the absolute value of potential energy 
as t~^, it is reasonable to expect that these quantities will become equal to each other at 
some instant of time, and next the plasma will evolve in the strongly nonideal regime. In 
such a case, let us consider a sufficiently small (but macroscopic) plasma volume where 
thermodynamic equilibrium is believed to be established and which, therefore, can be 
described by the multiparticle distribution function of the following general form: 



where r„ and v„ are the electron coordinates and velocities, R„ are the ion coordinates, 
and ^/ is the normalization factor. Since the kinetic energy of ions in the laboratory 

starting from the same initial temperature . Unfortunately, the plots in Fig. 3a of the above-mentioned 
paper are taken at the same "final" temperature (defined in some arbitrary instant of time). Most 
probably, this was done just to avoid merging the curves with a horizontal axis. As a result, it looks at 
the first glance that the theory considerably disagrees with the experiment at small rather than large 
times. On the other hand, it is written in the text of the article about the disagreement at large times, 
as should be expected from the physical formulation of the problem. 

^ Some theoretical arguments regarding suppression of the recombination in ultracold plasmas can be 

found, for example, in paper [12] and references therein. Anyway, even if the recombination is taken into 
account by the standard way [11], its influence on the general law of temperature evolution is quite small. 



2 Analytical Estimates 





n=l 
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experiments is much less than the electron kinetic energy, it will be ignored here; and, 
therefore, the electron motion will be treated at a fixed ion distribution. (If necessary, 
the same consideration can be easily performed when all kinds of the particles are taken 
into account.) 

Although the potential energy U in the regime of strong Coulomb's interaction has 
a very complex form and cannot be treated anymore as a small correction to the kinetic 
energy, the calculation of average value of some macroscopic quantity F depending only 
on the velocities v„ {e.g. kinetic energy) can be performed quite easily: 

/-(v)exp{-^[i:i^]}- 



(FM) = ' ' „ , (2) 



exp 
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dv 



because the integrals J exp{—U{r,'R.)/kBTe}dr in the numerator and denominator au- 
tomatically cancel each other. (For conciseness, v, r, and R denote here the sets of all 
velocities and coordinates of the electrons and ions, respectively.) 

Particularly, when the kinetic energy of an ensemble of particles is calculated, the 
integral in the numerator of formula ([2]) is reduced to the combination of Gaussian expo- 
nents, whose method of calculation is well known. By such a way, it can be shown that 
the average kinetic energy per one particle is given by exactly the same expression as for 
ideal gas: 

{k) = (3/2) kBT, ; (3) 

although, let us mention once again, it is valid for the plasma with arbitrarily strong 
Coulomb's interaction. 

Unfortunately, if we need to calculate the average values of quantities which are the 
functions of coordinates {e.g. the average potential energy), then the distribution func- 
tion ([T]) becomes actually useless, because the integrals involving the potential energy U 
in the exponent cannot be calculated in any reasonable approximation. Nevertheless, a 
special "way around" can be used here. Namely, let us relate the average potential energy 
per one particle (u) to the average kinetic one (k) by the virial theorem for the Coulomb's 
field [13], which is also valid at any strength of the interparticle interaction 



(k) = (1/2) |(«)| . (4) 



^ Let us mention that the virial theorem formulated for macroscopic bodies in some textbooks on 
statistical physics (see, for example, [14]) involves an additional term of the form 3PV, where P is the 
pressure, and V is the volume of the system. This term appears due to the surface integral for the 
particles interacting with a wall, where the potential energy is no longer the Euler homogeneous function. 
Since, in the case under consideration, the plasma is not confined by the walls, and the potential energy 
of interaction between its particles is everywhere the Euler homogeneous function, then no extra terms 
should appear in the virial theorem. 
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Of course, it is necessary to assume here the ergodicity of the system, i.e. that the 
quantities averaged over a statistical ensemble are equal to the ones averaged over time. 

Let us mention also that a necessary condition to apply the virial theorem is a 
bounded phase volume of the system, i.e. the particles must move in a limited spa- 
tial region with the velocities limited by the absolute value. The first condition, strictly 
speaking, is not satisfied for the plasma cloud infinitely expanding in space. Nevertheless, 
one can expect that the virial theorem can be reasonably applicable to the system whose 
phase volume is unbounded but increases with a small rate as compared to the charac- 
teristic velocities of its particles. This condition is well satisfied in the experiments with 
ultracold plasmas, because the characteristic time of variation in macroscopic parameters 
of the cloud (~ 10~^ s for the particular experimental setup [11, 15]) is much greater 
than the periods of microscopic motion of the electrons (10~^ -^ 10~^ s). We shall discuss 
the problem of applicability of the virial approximation in more detail in the end of the 
article. 

At the last step of our estimates, the average potential energy can be evidently 
expressed through the characteristic distance between the particles or the plasma density: 

(n) ~ e7(r) ~ e2r^^/^ (5) 

Finally, combining the formulas ©-([S]), we get TgOC n^^^. In particular, if the cloud 
expansion is inertial {i.e. linear in time) and, consequently, its concentration changes 
as t~^, then 

Teocr^ (6) 

Therefore, the presented model of "virialization" of the charged-particle energies in 
the regime of strong Coulomb's coupling well explains the both experimental features, 
namely: 

(a) the system "forgets" in the course of time about its initial temperature, i.e. the plasma 
clouds with various initial temperatures begin to evolve similarly; and 

(b) the particular form of the time dependence is close to instead of t~^, expected 
intuitively. The heat release due to inelastic processes (particularly, three-body recombi- 
nation) is not of importance here. 

3 Numerical Simulation 
3.1 Formulation of the Model 

To verify the above analytical estimates, we performed ab initio numerical simulation, 
based on the solution of the equations of classical mechanics for the multiparticle system. 
Our approach differs from the preceding works by the following items. 
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Figure 1: Sketch illustrating the interaction between an electron (e) and ion (i) both 
inside the basic cell (solid arrow) and with all its mirror images (dashed arrows). The 
ions are pictorially drawn by the large circles; and electrons, by the small ones (in fact, 
all the particles in our simulation were point-like) . 



The authors of the most of simulations of ultracold plasmas performed by now tried 
to include in their calculations as many particles as possible. Consequently, to keep the 
computational time in reasonable limits, they used a substantial simplification of the 
Coulomb's interactions, namely: 

(1) to avoid large integration errors during the close collisions, the Coulomb's potential 
was cut off of smoothed out at the small distances {e.g. 1/r was changed to l/(r + ro) or 
something like that) [12, 16, 17]; 

(2) to improve convergence of the Coulomb sums at large distances, the researchers used 
the multipole-tree method [18] and various approximations for the equations of motion of 
the light particles (electrons), such as the particles-in-cell (PIC) method [19] or "Vlasov 
approximation" for the electron component of plasma [20] {i.e., actually, the introduction 
of self-consistent field, ignoring the interparticle correlations). 

All these approximations distort the functional dependence of the Coulomb's poten- 
tial so strongly that one can hardly expect the virial relations to be satisfied, because the 
virial theorem is very sensitive to the particular form of the potential energy. 

As distinct from the above approaches, we tried in our simulation to integrate the 
equations of motion of the charged particles in the real electric microfields as accurately 
as possible, without any artificial distortion. With this aim in view, we used the basic 
cell with a relatively small number of particles {e.g., a few dozens) which was assumed 
to be supplemented in all directions by infinite number of mirror cells. Coulomb forces 
acting on each particle in the basic cell were calculated taking into account not only the 
particles in the same cell but also in all its mirror images, until the specified accuracy is 
achieved (Fig. 1). 
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So, we worked with a relatively small number of the equations of motion — only for 
particles in the basic cell. As a result, we could perform integration with a very small step 
and, therefore, automatically avoid the problem of close collisions without any artificial 
modification of the Coulomb's potential. On the other hand, the total number of the 
charged particles participating in Coulomb interactions (including the mirror images) was 
in our calculations between 100 000 and 3 000 000. This is not less and even usually much 
greater than in the simulations by other researchers. 

3.2 Initial Equations 

As was already mentioned in Sec. 121 the ion kinetic energy in the experiments with 
magneto-optical traps is usually small as compared to the kinetic energy of electrons. 
Consequently, the ions in our simulation were assumed to move from the very beginning 
by the inertial {i.e., linear in time) law; while the electron dynamics was described taking 
into account Coulomb's interactions not only inside the basic cell but also with infinite 
number of its mirror images: 



Here, Rj {i = 1, A^) are the ion coordinates, {i = 1, ZN) are the electron coordi- 
nates, A^ is the number of ions in the basic cell, Z is the ion charge (in the calculations 
presented below in this article, the ions were singly charged, i.e. Z = 1), m and e are 
the mass and absolute value of the electron charge, L is the linear size of the basic 
cell, (Z = 1,2,3) are the unit vectors of the Cartesian coordinate system, and sub- 
scripts Hk {k = 1,2,3) specify the number of the mirror cell in each direction. Let us 
mention that a software code actually performed summation over the mirror cells not 
from — oo to +oo, as in formula ([7]), but starting from the central (basic) cell, one shell 
of cells after another, until a specified criterion of convergence of the Coulomb sums is 
satisfied. 

Yet another well-known problem in modeling the freely-expanding plasmas is a con- 
siderable variation in the spatial scale of the system (and, consequently, in the values of 
Coulomb's forces) during its evolution. As a result, it becomes difficult to choose the 
method of numerical integration ensuring a stable accuracy for the entire solution. In our 
modeling, this problem was resolved by introduction of a "scalable" coordinate system. 






(7) 
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expanding in space with a mean plasma expansion velocity. In other words, size of the 
basic cell was taken to be increasing linearly in time: 



Lo + uot (8) 

(which corresponds, from the physical point of view, just to the incrtial motion), and the 
coordinates of all particles were normalized to this timc-dcpcndcnt scale. 

The ions, moving by inertia, will be exactly at rest in such coordinates, while the 
electrons will move within a cell of the fixed size. When the electron crosses one of 
boundaries of the cell, it is assumed to appear at the opposite side, as it is usually done 
in the method of molecular dynamics. 

At last, to get the final equations of motion for the electrons, which will be used in 
the simulation, all physical quantities should be reduced to the dimensionless form. Let 
the time- dependent unit of length / be the characteristic interparticle distance, determined 
by the following way. The total number of particles in the basic cell equals (Z + 1)A^; 
so the volume per one particle will be L^/{{Z + 1)N). Consequently, the characteristic 
linear size can be determined as 

(Z+l)V3Ari/3 ■ (9) 

Particularly, at the initial instant of time (which from here on will be denoted by sub- 
script 0) we get: 

^0 = (Z + l)V3iVV3 ■ (10) 
The characteristic time scale r can be formally introduced, for example, by the virial 
relation taken at the initial instant of time: (1/2) m (Iq/t) ^ = (1/2) Ze^/ Iq , from which 
we obtain 

"^^'''/o^/l (11) 



/ m \ 



Within a numerical factor, this coincides with the inverse Langmuir frequency (as could 
be expected from the dimensionality arguments). 

Prom here on, the quantities normalized to I and r will be denoted by asterisks. It 
can be easily shown that the physical scale I will be a function of the dimensionless time 
of the form: 

h{l+ult*), (12) 

where 

uI = Uot/Lq; (13) 
while a dimensionless length of the basic cell is constant: 

L* = {Z + iy/^N^/\ (14) 
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Therefore, in the variables specified above, equations of the electron motion ([7j) will 
take the form: 

r:+2u;{l+ultT'i:= (l+<0-'F*, (15) 

where dot denotes a derivative with respect to the dimensionless time t*, and F* is the 
total Coulomb's force acting on ith electron from all other electrons and ions: 

^ "R* y* 1 „* v.* 

* Z^|"R*_r*|3 7 Ir*— r*|3 

As follows from equation (|T5l) . the effect of inertial plasma expansion in the "ex- 
panding" coordinate system looks like the influence of an effective dissipative force, which 
is proportional to the electron velocities. Therefore, temperature of the electron gas re- 
sults from the balance of two effects — on the one hand, acceleration and heating of the 
electrons due to Coulomb's interactions and, on the other hand, their deceleration and 
cooling by the above-mentioned dissipative forces. 

3.3 Method of Computation 

To write conveniently the subsequent formulas, let us introduce the auxiliary quantity: 

s(r) = (1 + M*r)-^ (17) 

Then, following the standard procedure, the second-order equation (ITB]) can be 
rewritten as a set of two equations of the flrst order: 

rr=v*, (18a) 
v* = -2<sv*+ s^F*, (18b) 

which can be solved by any available method of numerical integration. We used Runge- 
Kutta method of the second order. 

The set of equations fll8ap and (]18bp is solved in the region 



- L72 < (r*)i < , where i = 1, ZN and / = 1, 2, 3 (19) 

(the subscript / denotes the number of the Cartesian coordinate) with standard molecular- 
dynamic boundary conditions: a particle leaving the basic cell of simulation through one 



9 



of its sides is replaced by the particle entering the cell with the same velocity through the 
opposite side. The initial coordinates of electrons are given by the random-number gen- 
erator as a uniform statistical distribution over the region flTI?]) : and the initial velocities, 
as Gaussian distribution with a specified dispersion. 

As was already mentioned above, an ion motion in the simulated experimental con- 
ditions can be considered as inertial; so that the normalized ion coordinates in the "ex- 
panding" reference system remain constant. At the initial instant of time, they are taken 
as a uniform statistical distribution over the region 

- L*/2 < {R*)i < L*/2 , where i = l,...,N and 1 = 1,2,3. (20) 

It is unnecessary, of course, to specify the initial velocities and to solve the equations of 
motion for the ions. 

Finally, let us discuss in more detail how to calculate a kinetic energy of the electron 
motion K with respect to the mean plasma flow, which is related to their temperature by 
formula ([3]). We can see here one more advantage of the introduced coordinate system: the 
required kinetic energy is calculated in this system just by differentiating the normalized 
(dimensionless) electron coordinates with respect to time; and one should not differentiate 
the normalization factor l{t) itself, because its variation is associated with a kinetic energy 
of the plasma motion as a whole: 

--TE(^.)>i^E(J.^;J-i^E(v.-)^ (21) 

i=l V / rcl j=i \ / j=i 

By introducing the normalization factor 

^■-j'i- (22) 

we find that the dimensionless kinetic energy of the relative motion is given by the ex- 
pression: 

ZN 

= s'' (^*) ' • (23) 



i=l 

Just this formula (after division by the total number of particles in the basic cell) will be 
used below to calculate the plasma temperature. 



3.4 Particular Case of the Ideal Gas 

Before presentation of the results of numerical modeling, let us discuss one particular 
example, which can be completely solved in analytic form and demonstrates a self- 
consistency of the approach used. Namely, let us consider the case of an ideal gas, when 
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the forces of interparticle interaction F* disappear at all (for example, because the electric 
charges tend to zero). Then, equation ( 118b|) is reduced to 

vlk = -2 uls v*, (24) 

(subscript i denotes here the number of particle, and k is the number of coordinate). 
After the integration of this differential equation, taking into account the particular form 
of function s{t*) given by formula f|T71) . we get: 

vlk = «.)o (1 + <tT' = «fc)o s'in ■ (25) 

At last, substituting solution ( 125|1 to general expression for the kinetic energy (123|1 . 
we find: 

ZN 

^''E«)o^'= ^0^'(O- (26) 

i=l 

Therefore, the electron temperature oc K cc K* will change in time by the law: 

Te{t) oc s^(r) oc (1 + M*r)~2 cx t'^ at large t, (27) 
as should be expected for the inertial expansion of the ideal monatomic gas. 

3.5 Basic Parameters of the Numerical Simulation 

Before presentation of the results of numerical simulation for strongly-coupled plasmas, 
let us describe in more detail the basic parameters used in the computation presented 
below in Figs. 2 and 3. 

The number of particles in the basic cell: A^tot = (2^ + 1)A^ = 20 {i.e., 
10 electrons and 10 ions). Due to the so small number of particles, we were able to 
perform integration with high accuracy (see below for more details) over the time interval 
Atot^* ~ 1000. The total number of particles taken into account in the calculation of 
Coulomb sums (including the mirror cells) was, on the average, a few hundred thousand 
(more exactly, 137180 to 3 327500). 

In other versions of simulation, for example, when the number of particles in the 
basic cell was increased by an order of magnitude — up to Atot = 200, we were enforced to 
decrease the interval of integration by 3-^4 times; but the behavior of macroscopic plasma 
parameters remained qualitatively the same. When Atot was further increased by a few 
more times, the reachable interval of integration became so short that it was difficult to 
draw reliable conclusions on the laws of evolution of the plasma parameters with time. 

The step of integration and the accuracy of calculation of the Coulomb's 
forces: At* = 10~'^-^ 10"^ and ep = 10~^-i- 10^^. These values were chosen empirically, 
by a series of test simulations. Starting from the above-written values, the behavior of 
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macroscopic plasma parameters became reproducible, i.e. independent of further increase 
in the accuracy of computation. Let us mention that the integration step At* = 10^^ was 
usually sufficient. However, at some "unfavorable" initial distributions of theparticles, 
there might be a few close collisions when such integration step was too largejj In these 
cases, it should be reduced to 10"^ and sometimes even smaller. 

Velocity of the inertial plasma expansion: Uq = 0.1 . Such value of the expan- 
sion velocity at the time interval under simulation results in increasing the plasma cloud 
size by a few dozen times, which corresponds to the experimental conditions [11]. 

Root-mean-square scatter of the initial velocities of the particles: cr*o = 3.0 
(in each Cartesian coordinate), i.e. initially the plasma was taken to be slightly nonideal. 
In some other simulations, we used greater values of cr*Q, i.e. started from a more ideal 
plasma state. In such cases, the initial stage of evolution of the electron temperature was 
close to the expected dependence Te oc t~^. In subsequent modeling, to avoid spending 
a lot of time for the integration over the physically trivial interval, we preferred to start 
just from CT*Q = 3.0 . 

3.6 Results of the Simulation 

The results of our numerical simulation for the electron temperature are presented by 
crosses in logarithmic scale in Fig. 2. After excluding from a consideration the earliest 
period of time, when relaxation processes occurred, the points in the remaining time 
interval lie almost along a straight line, corresponding to the power law Tg oc t". The 
respective exponent was found to be in the range a = —(1.08-^1.25). This is quite close to 
the value a = —1, following from the simple virial estimate, and is in perfect agreement 
with the experimental value a = — (l.l-f-1.3) [11]. 

It is reasonable to ask if the dependence obtained is really caused by the effect 
of virialization, discussed in Sec. [2p To answer this question, we plotted in Fig. 3 the 
temporal dependences of average kinetic energy and one-half the absolute value of the 



potential (Coulomb's) energj.^! in ordinary (not logarithmic) coordinates. These quantities 
were obtained by the method of moving average over the interval Aavr^* = 20 {i.e. 10 units 
of dimensionless time in both sides from the center). It is seen in this figure that at large 
time, t* ~ 500, the curves k*{t*) and \u*{t*)\/2 begin to coincide with a sufficiently good 
accuracy, about 10%. On the other hand, it is important to emphasize that the above- 



^ The insufficiently small value of the integration step in close collisions manifests itself, first of all, 
as a "sling effect", i.e. a sudden ejection of the electron during a passage of pericenter of the orbit. In 
macroscopic plasma description, such situations look as non-physical saw-toothed peaks in the dependence 
of kinetic energy on time. 

^ Let us mention that the value of Coulomb's energy used here is not just an estimate by the char- 
acteristic interparticle distance, but it is the quantity accurately calculated by the summation over all 
particles, taking into account the mirror cells, as in the case of Coulomb's forces. 
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\nt* 

Figure 2: Electron temperature obtained by the numerical simulation as function of time 
in the logarithmic coordinate system. The inclined straight lines show the power-like 
dependences of the form Tg oc t"; the values of a being presented near the respective 
lines. 

mentioned power law of evolution of Te with exponent a = — (l.l-=-1.3) is established 
much earlier, when the virialization criterion is satisfied only on the order of magnitude. 
Probably, this is just the reason why a differs from the exact value of —1. This question 
requires a further careful study. 

At last, it is interesting to compare our findings with the results of paper [21], which 
was based on a quite complex "hybrid model" of ultracold plasmas. It involved the quasi- 
hydrodynamic equations for ions, description of the electron dynamics by Monte-Carlo 
method, as well as additional inclusion of some collisional processes by the a priori in- 
formation about their cross-sections. As a result, it was in particular found that, if the 
three-body recombination was not included explicitly, then the electron coupling parame- 
ter Fe increases infinitely with time at any initial temperature. On the other hand, if the 
effects of three-body recombination and subsequent inelastic collisions of electrons with 
the produced Rydberg atoms were included to the equations as additional terms, then the 
coupling parameter was stabilized in the course of time at the level Fg ~ 0.2 (see Fig. 3 
in paper [21]). 
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Figure 3: Temporal behavior of the average kinetic energy per one particle k* (curve 1) 
and one-half the absolute value of average potential (Coulomb's) energy \u*\/2 (curve 2). 



Unfortunately, this conclusion does not agree with the results of paper [17], pub- 
lished almost at the same time. Its authors used a "pure" molecular-dynamic model, 
without any explicit inclusion of the three-body recombination, and found a stabilization 
of Fe at the level about unity. However, it should be mentioned that only the slightly- 
expanding plasmas were considered. As regards our simulation, it also does not take 
into consideration any special corrections for the recombination and gives the asymptotic 
values of coupling parameter Fg ~ 1, i.e. in agreement with [17], also in the case of very 
large plasma expansion (by a few dozen times, in terms of the linear size)^ 

Let us emphasize that, if the recombination processes are not taken into account 
explicitly in the equations of molecular dynamics, this does not mean that these processes 
are completely ignored. In fact, our model already involves their description from the 
first principles: as follows from a more careful analysis of the computational results, 
the scattering of crosses in Fig. 2 and spiky-like behavior of the curves in Fig. 3 are 
caused just by formation of quasi-bound states of the electrons and ions {i.e. by the 
onset of recombination between the charged particles). If the electron-ion pair with a 
sufficiently large eccentricity of the orbit is formed, then every passage of the electron in 
the vicinity of the ion will be associated with a sharp outburst of both kinetic energy and 
the absolute value of potential energy. This looks as a series of sharp equidistant peaks in 
macroscopic parameters of the plasma. A few series of such peaks were clearly observed in 



^ Comparing the values of coupling parameter T obtained in the various works, one should keep in 
mind that definitions of this parameter may be slightly different, by a factor about unity. 
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our simulations. However, consideration of this phenomenon requires a separate article, 
and it will not be discussed here in more detail. 

4 Conclusions 

1. The law of evolution of the electron temperature close to can be explained, in 
the first approximation, by a simple analytical model based on the virialization of energy 
of the charged particles. 

2. The numerical simulation from the first principles, using the scalable coordinate 
system and accurately taking into account Coulomb's interactions, results in some devi- 
ation of the exponent, namely, from —1 to — (1.08-i-1.25), which is in perfect agreement 
with the experimental value. 

3. Therefore, the decelerated law of decrease in the electron temperature is, first of 
all, a manifestation of the specific equation of state of the cold nonideal plasma rather 
than a result of additional heat release due to recombination. 

4. It was unexpectedly found in the numerical simulation that the law of variation in 
the electron temperature of the form Tg oc i-(i 08-^i-25) jg estabhshed very quickly — when 
the virial relation for energies is satisfied only within a factor on the order of unity. 

5. A more careful analysis of the results of our simulation reveals also some finer 
effects, for example, the quasi-periodic oscillations of energy caused by the formation of 
bound electron-ion states. 
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